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Programming techniques which extend the capabilities of FeynArts and FormCalc are introduced and explained 
using examples from real applications. 



1. Introduction 

FeynArts [T] and FormCalc 0j are programs 
for the generation and calculation of Fcynman 
diagrams. The main problem they solve is the 
computation of the cross-section for a given pro- 
cess, and indeed this is meanwhile a heavily au- 
tomated procedure that requires little user input. 
Feynman diagrams have a lot of other uses, how- 
ever, and this note shows some examples of how 
to extend the capabilities of FeynArts and Form- 
Calc to produce results not within the 'standard 
canon. ' 

Although cash-strapped institutes may tend to 
believe otherwise, it is actually a feature, not 
a bug, that FeynArts and FormCalc are largely 
written in Mathematica, for it is precisely the 
availability of a powerful language that makes it 
easy for the user to examine and modify results, 
and also to extend the functionality of these sys- 
tems. Mathematica is particularly suited for the 
latter as it is technically an expert system, where 
knowledge is added in the form of transformation 
rules. Thus one usually does not have to modify 
the existing program code (possibly including re- 
compilation), but just add new rules. Sections 
[3J 01 and [3] show some real-life examples of such 
extensions. Section H3 exhorts on the concept of 
abbreviations in FormCalc, where some impor- 
tant improvements have been made recently. 

2. Programming Diagram Filters 

Question: What if FeynArts' diagram selection 
functions are not enough? What if I want, say, 



only diagrams with a fermion loop? 

Answer: Inspect the internal structure of the 
inserted topologies, i.e. the InsertFields out- 
put. The outermost structure is a TopologyList: 

TopologyList [ ] Ui, t 2 , ...] 

Contained in this are topologies U of the form 

Topology [_] [ ] -> 

Insertions [Generic] [pi, gi, ...] 

The generic insertions gi specify the fields running 
on each line of the topology, e.g. 

Graph [__] [Field [1] -> F, ...] -> a 

Now recall that FeynArts distinguishes three lev- 
els of fields (generic, classes, particles) where the 
field's space-time properties are fixed at the low- 
est (generic) level [3]. 

The fermion-ness of a field is a space-time prop- 
erty, so we need not go deeper than the generic 
level. In order to look at generic-level diagrams, 
the selection function for fermion loops must first 
"propagate down" into the generic level: 

FermionLoop [ t : TopologyList [ ] [ ] ] : = 

FermionLoop/® t 

FermionLoop [ (top : Topology [_] [ ] ) -> 

ins : Insertions [Generic] [ ] ] : = 

top -> TestLoops [top] /<§ ins 

The subsidiary function TestLoops selects only 
the fields carried by propagators in the loop: 

TestLoops [top_] := TestFerm[ Cases [top, 
Propagator [_Loop] [_ , _, field_] -> 
field] ] 
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and hands on to TestFerm which tests if, for a 
concrete diagram, the selected fields are indeed 
all fermions: 

TestFerm [f ields_] [gi_ -> ci_] := 
(gi -> ci) /; 

MatchQ [fields /. List®® gi , {F..}] 

TestFerm [_] [_] := Sequence [] 

Note this typical (for Mathematica) construct: 
the first instance handles diagrams positively 
identified as having a fermion loop, and the sec- 
ond one is the 'fall- through' for all other cases. 

3. Tweaking Model Files 

Question: Is there a good way to make (small) 
changes to an existing model file? 

Answer: It is bad practice to copy the model 
file and modify the copy, for two reasons: 

1 . It is typically not very transparent what has 
changed. 

2. If the original model file changes (e.g. due to 
bug fixes), these do not automatically prop- 
agate into the derivative model file. 

Better: Create a new model file which reads the 
old one and modifies the particles and coupling 
tables. To this end one needs to know that the 
model file defines two main objects: 

• M$ClassesDescription is the list of parti- 
cle definitions, 

• M$CouplingMatrices is the list of couplings 
of the form C [fields] == expr. 

As an example consider introducing enhance- 
ment factors for the b-b-ho and b-b-Ho Yukawa 
couplings in the MSSM. The new model file is 
both compact and makes it immediately obvious 
what has changed: 

Block [{$Path = $ModelPath}, « MSSM. mod] 

EnhCoup [ (lhs : C [F [4, {g_ , _}] , -F [4 , _] , 

S[h:l|2]]) == rhs_] := 
lhs == Hff[h,g] rhs 
EnhCoup [other_] = other 

M$CouplingMatrices = 

EnhCoup/® M$CouplingMatrices 



Note that the enhancement factors depend on 
the fermion generation. This is because the cou- 
plings are defined at the classes level, i.e. for the 
class F[4] of down- type quarks {d, s,b}, not for 
the bottom quark alone. Thus one needs to set 
Hff [h,l] = Hff [h,2] = 1. 

A printout of the new Feynman rules can be 
obtained with the WriteTeXFile .m program that 
comes with FeynArts. 

4. Scripting Mathematica 

Question: How can I do efficient batch process- 
ing with Mathematica? 

Answer: Put everything into a script, using sh's 
Here documents: 

#! /bin/sh Shell Magic 

math « \_E0F_ start Here document 

« FeynArts' 

« FormCalc' 

top = CreateTopologies [. . .] ; 

_E0F_ end Here document 

Everything between "« \tag" and "tag" goes to 
Mathematica as if it were typed from the key- 
board. Note the "\" before tag, it makes the shell 
pass everything literally to Mathematica, with- 
out shell substitutions. This is important because 
Mathematica uses many characters which have a 
special meaning to the shell, such as $, [, {, etc. 
The advantages of this method are: 

1. Everything is contained in one compact 
shell script (i.e., a text file), even if it in- 
volves several Mathematica sessions. 

2. Such a script can easily be run in the back- 
ground, or combined with utilities such as 
make. 

3. One can seamlessly combine Mathematica 
and shell programming. A slightly subtle 
issue is how to get shell variables such as 
command-line arguments into Mathematica 
even though substitutions in the Here file 
are turned off through the use of \ in \tag. 
Solution: pass them on the command-line 
when invoking Mathematica: 
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#! /bin/sh 

math -run "argl=$l" . . . « \END 
END 

Debugging hint: the -x flag makes sh echo every 
statement, it can be added after the Shell Magic: 

#! /bin/sh -x 

5. Not the Cross-Section 

Question: Can I get things out of FormCalc 
other than the cross-section? Can I, for example, 
compute the Wilson coefficients for b — > 57? 

Answer: The relevant operators for b — > S7, the 
prefactors of which are the Wilson coefficients, are 

O7 = J^™b (sj\oJ+a^ \bi) SijF^ (1) 

° 8 = lf^ mb ^ T % G « (2) 

Generating the partonic diagrams and amplitudes 
is of course a standard exercise: 

tops = CreateTopologies [1 , 1 -> 2] 
ins = InsertFields [tops , 

F[4,{3}] -> {F[4,{2}] ,V[1]>] 
amp = CalcFeynAmp[CreateFeynAmp[ins] , 

FermionChains -> Chiral] 

For the 0$ coefficient, the photon V[l] has to be 
replaced by the gluon V[5] . The FermionChains 
option ensures that the result contains the chiral 
Dirac chains needed to read off the Wilson coef- 
ficients, rather than the default Weyl chains. 

CalcFeynAmp collects Dirac chains and colour 
matrices in a function Mat. For identifying the 
operators it is thus simplest to construct a func- 
tion id that replaces Mat. 

The identification of the operators is slightly 
nontrivial because FormCalc strives to reduce the 
number of Lorentz indices as much as possible 
and thus turns uj±^^ v into u!±j^. This can be 
reversed by applying the Gordon identity 

/I u I \ (Pi +P2T , , 1 , . 

(s 2 \uj±T \ s i) = ( s 2 w± Si) + 

mi 

i(P2 -Pl)u , 1 in, 1 \ , Q . 

{s 2 \ u) ± a^ v |si) . (3) 



Since this is contracted with e^, the photon/gluon 
polarization vector, we can use the momentum- 
space correspondence a^F^ — > 2a^(p 2 — Pi)^e^ 
to identify the F^-term, thus 

2 

(s 2 | \si) = — e-pi (s 2 | u)± \si) + 
mi 

-2- {s 2 \uj±ct» v \ Sl ) F^ . (4) 

id[r_. DiracChain[s2_Spinor, om_, eps_, 
sl:Spinor [pl_, ml_, _] ] ] := 
2/ml r Pair[eps, pi] * 

DiracChain[s2, om, si] + 
1/(2 ml) id[r sig[om]] 

It is not even necessary to explicitly write out the 
cr^-tcrm since its only purpose is to match the 
operators: 

id[r_. sig[om_] SUNT[i_, := 
r 07[om]/(EL MB/ (16 Pi~2)) 

id[r_. sig[om_] SUNT[a_, i_, j_]] := 
r 08[om]/(GS MB/ (16 Pi~2)) 

Now we can apply this function to the amplitude: 

amp = Plus@@ amp //. Abbr [] /. Mat -> id 

c7 = Coefficient [amp, 07 [6]] 
c8 = Coefficient [amp, 08 [6]] 

Using FormCalc's output functions it is also 
pretty straightforward to turn these expressions 
into Fortran code: 

file = OpenFortran ["bsgamma. F"] 

WriteString [file , 

SubroutineDecl["bsgamma(C7, C8)"] <> 
"\tdouble complex C7, C8\n" <> 
"#include V'model .h\"\n" <> 
" #include \ " loopt ools . h\ " \n" ] 

WriteExpr[f ile, {C7 -> c7, C8 -> c8}] 

WriteString [file, "\tend\n"] 

Close [file] 

6. Abbreviations 

The automated introduction of abbreviations is 
one of the key concepts in FormCalc. It is crucial 
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in rendering an amplitude as compact as possible. 
The main effect comes from three layers of recur- 
sively defined abbreviations, introduced when the 
amplitude is read back from FORM, i.e. during 
CalcFeynAmp. For example: 



10-30% speed-up. It works in two steps. First, 
redundant parts are removed, e.g. the abbrevia- 
tions 



AbbSum29 = Abb2 + 



Abb22 



+ Abb23 + Abb3 



Abb22 = Pairl Pair3 Pair6 



|Pair3 = Pair [e [3] , k[l]] 

Written out, this abbreviation is equivalent to 

Pair[e[l] , e[2]] Pair [e [3] , k [1] ] Pair [e [4] , k [1] ] + 

Pair [e [1] , e [2] ] Pair [e [3] , k [2] ] Pair [e [4] , k [1] ] + 

Pair [e [1] , e [2] ] Pair [e [3] , k [1] ] Pair [e [4] , k [2] ] + 

Pair[e[l] , e[2]] Pair [e [3] , k [2] ] Pair [e [4] , k[2] ] 

In addition to these abbreviations assigned by 
CalcFeynAmp, FormCalc introduces another set of 
abbreviations for the loop integrals when gener- 
ating Fortran code, i.e. during WriteSquaredME. 

6.1. Categories 

Both of the aforementioned types of abbrevi- 
ations, but in particular the latter, are costly in 
CPU time. It is thus key to a decent performance 
that the abbreviations are grouped into different 
categories: 

1. Abbreviations that depend on the helicities. 

2. Abbreviations that depend on angular vari- 
ables. 

3. Abbreviations that depend only on ^/s. 

Correct execution of the different categories guar- 
antees that almost no redundant evaluations are 
made. For example, for a 2 — > 2 process with 
external unpolarized fermions, statements in the 
innermost loop over the helicities are executed 2 4 
times as often as those in the loop over the an- 
gle. This technique of moving invariant expres- 
sions out of the loop is known as 'hoisting' in 
computer science. 

6.2. Common Subexpressions 

Another optimization method, common subex- 
pression elimination, is implemented in the func- 
tion OptimizeAbbr and can yield some additional 



- Abb74 - Abb87 
+ Abb74 - Abb87 



AbbSum2 -> Abb9 - Abb87 

AbbSum8 -> Abb9 - Abbl3 

AbbSum9 -> Abb9 + Abbl3 

are replaced by 



AbbSum2 -> Abb9 - Abb87 

AbbSum8 -> AbbSum2 - Abbl3 - Abb74 

AbbSum9 -> AbbSum2 + Abbl3 + Abb74 

In the second step, common parts are put into 
temporary variables, thereby simplifying the last 
lines further to 

AbbSum2 -> Abb9 - Abb87 
helpl -> Abbl3 + Abb74 
AbbSum8 -> AbbSum2 - helpl 
AbbSum9 -> AbbSum2 + helpl 

Both optimization techniques together make the 
generated code essentially as fast as hand-tuned 
code. 

6.3. The Abbreviate Function 

The new Abbreviate function extends the ad- 
vantages of the abbreviation system to arbitrary 
expressions. Its usage is for example: 

abbrexpr = Abbreviate [expr , 5] 

The second argument, 5, determines the level be- 
low which abbreviations are introduced. The level 
determines how much of expression is 'abbrevi- 
ated away' i.e. how much of the structure is pre- 
served. In the extreme, for a level of 1, the result 
is just a single symbol. Abbreviationing also has 
the 'side effect' that duplicate expressions are re- 
placed by the same symbol. 

This new type of abbreviations for subexpres- 
sions has to be retrieved separately from the other 
ones with Subexpr [] . 

The most important option of Abbreviate is 
Preprocess. It is used e.g. as follows: 

abbrexpr = Abbreviate [expr , 5, 
Preprocess -> Simplify] 
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and specifies a function which is applied to each 
subexpression before the abbreviations are intro- 
duced. 

At some 30 sec. execution time for Abbreviate, 
the typical speed-up was a factor 3 in MSSM cal- 
culations. 

6.4. Auxiliary functions 

When simplifying large expressions, it is of- 
ten desirable to have a replacement for Simplify 
which is faster by performing only specific sim- 
plifications. This is because Simplify is quite ef- 
ficient on short expressions but increasingly slow 
on longer ones. 

DnSize constructs a special function for simpli- 
fication that does different things depending on 
the size of its argument. For example, 

f = 0nSize[100, Simplify, 

500, DenCollect, 
Pool] 

If the Leaf Count of the argument of f is 

• below 100, Simplify is used, 

• between 100 and 500, DenCollect is used, 

• above 500, Pool is used. 

DenCollect collects terms with denominators 
that are identical up to a numerical constant. 

Pool combines terms with common factors. 
Unlike Factor, it looks at the terms pairwise and 
can thus do ab + ac + d — » a(b + c) + d fast. Pool 
will not factor out very small expressions because 
the effect of this on the size of the whole expres- 
sion is typically not worth the effort. 

7. Summary 

We presented several examples of programming 
techniques which enable the user to extend the 
capabilities of FeynArts and FormCalc with com- 
paratively little effort: 

• programming diagram filters, 

• applying (small) changes to model files, 

• scripting Mathematica, 

• computing quantities other than the cross- 
section. 



In addition, we elucidated the important concept 
of abbreviations in FormCalc which is extended 
to arbitrary expressions by the new Abbreviate 
function. 
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